//      gene_processing.h
//
//      Basic sequence processing, read/write genome, annotation file
//      
//      Copyright 2012 tqlong <tqlong@espada>
//      
//      This program is free software; you can redistribute it and/or modify
//      it under the terms of the GNU General Public License as published by
//      the Free Software Foundation; either version 2 of the License, or
//      (at your option) any later version.
//      
//      This program is distributed in the hope that it will be useful,
//      but WITHOUT ANY WARRANTY; without even the implied warranty of
//      MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
//      GNU General Public License for more details.
//      
//      You should have received a copy of the GNU General Public License
//      along with this program; if not, write to the Free Software
//      Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston,
//      MA 02110-1301, USA.

#ifndef GENE_PROCESSING_H
#define GENE_PROCESSING_H

#include <string>
#include <vector>
#include <boost/program_options.hpp>
#include "stl-helper.h"

/** A genome which can be very long string of ACGT nucleotides */
class Genome : public std::string
{
  typedef std::string __BASE;
 public:
  Genome();
  Genome(const std::string& s);

  void read_genome(const std::string& genomeFile);
 protected:
  std::string firstLine;
};

/** An annotation entry of a gene (including upstream sequence) */
struct _annotation 
{
  std::string gene_id;
  std::string gene_name;
  int         gene_number;
  std::string gene_direction;
  int         gene_start;
  int         gene_stop;

  std::string gene;
  std::string gene_upstream;

  void print();
  int length() const { return gene_stop - gene_start + 1; }
};

/** A list of sequences that has been extracted from the genome 
  * This is the main data class for gene algorithms
  */
class Sequences : public std::vector<std::string>
{
  typedef std::vector<std::string> __BASE;
 public:
  void count_occurence(const std::string& pattern, std::vector<int>& count) const;
  int max_length() const;
 public:
  static void check_occurence(const std::string& pattern, const std::string& seq, std::vector<int>& pos);
  static bool compare(const std::string& pattern, const std::string& seq);

  void read_sequence(const std::string& seqFile);
 public:
  std::vector<std::string> annotation_;
};

/** Sequences that are divided into groups */
class GroupSequences : public std::vector<Sequences>
{
public:
  void group_sequences(const Sequences& seqs);
  GroupSequences removeSharp() const;
};

/** Read annotation file functionality 
  * Write upstream sequences
  */
class Annotation : public std::vector<_annotation>
{
  typedef std::vector<_annotation> __BASE;
 public:
  void read_annotation(const std::string& annotateFile);
  void read_annotation_simple(const std::string& annotateFile);
  void read_annotation_simple1(const std::string& annotateFile);
  void get_gene_sequence(const Genome& genome, int n_upstream = 50);

  void save_upstream(const std::string& upstreamFile);
  void all_upstream(Sequences& seqs);
 protected:
};

class PatternGenerator
{
 public:
  PatternGenerator(boost::program_options::variables_map& vm, const Sequences& seqs, const std::string& entropyFile);

  void run();

  static std::vector<std::string> next_candidate(const std::string& pattern);
 protected:
  boost::program_options::variables_map& vm_;
  const Sequences& seqs_;
  std::string entropyFile_;
  std::vector<std::string> pattern_;
};

class MotifGibbsSampler
{
 public:
  MotifGibbsSampler(boost::program_options::variables_map& vm, const Sequences& seqs);

  void run();
  const vi_t& pos() const { return motif_position_; }
  const vvd_t& freq() const { return freq_; }
  void compute_likelihood(int seq_idx, vd_t& likelihood);
  double log_likelihood();

 protected:
  void prepare_count();
  void gibbs_sampler(int seq_idx);
  void count(int seq_idx);
  void compute_freq(int seq_idx);
  void sample(int seq_idx);
  void print_debug(int iter);
  void phase_shift();
  void shift(const vi_t& current_pos, vi_t& pos, int dpos);
  void compute_max_position();

  boost::program_options::variables_map& vm_;
  const Sequences& seqs_;
  std::string acgt;
  std::string outputFile_;
  std::string freqFile_;
  std::string debugFile_;
  int max_iter_;
  int motif_length_;
  vi_t motif_position_;
  int prev_idx_;

  vvi_t count_;
  vvd_t freq_; 
  std::vector<vi_t> pre_count_;
};

class MotifStats
{
 public:
  MotifStats(boost::program_options::variables_map& vm, const Sequences& seqs);
  void run();
  void output();
 protected:
  void count_stats(const std::string& spacer);
  void reset_stats();

  boost::program_options::variables_map& vm_;
  const Sequences& seqs_;
  std::string acgt;
  std::string motif_;
  std::string outputFile_;

  vd_t background_freq_; // ACGT frequency in spacer
  vvd_t base_pos_freq_;  // ACGT position frequency in spacer
  vvd_t pos_base_freq_;  // ACGT frequency in each position in spacer
  vd_t length_freq_;     // spacer length frequency
  double count_;
};

#endif
